Consuming GIS REST Services in R
Section Contact:Rex Robichaux (rex.robichaux@deq.virginia.gov)
GIS Rest Service Data
DEQ Makes a wide variety of it’s geospatial data available in the form of published REST services (Esri Map Services). First- a little background information on GIS Services in general:
Esri ArcGIS REST API documentation
One advantage of learning to use Esri REST services over these other formats is that a user may filter their data prior to storing it locally, limiting download size and required data cleaning/filtering. Users may also go on to implement use of these APIs within programs and scripts and know that they are using the most current data available. Another advantage of utilizing Esri REST platforms is that they are available to non Esri clients allowing users to avoid costly subscription fees surrounding ArcGIS and other Esri software.
Primary DEQ REST Endpoints:
Internal (Staff Only) Services: https://gis.deq.virginia.gov/arcgis/rest/services/staff
External (Public Facing) Services: https://gis.deq.virginia.gov/arcgis/rest/services/public
It’s important to note that all GIS Services have an API access point to perform queries, which can greatly assist in scripting/query development.
Once you have located a service, and a particular layer within a service (denoted by a /# format (ending in a number) such as “https://gis.deq.virginia.gov/arcgis/rest/services/staff/DEQInternalDataViewer/MapServer/72” for TMDL Watersheds), you can find a “Query” operation at the bottom of the page under the “Supported Operations” section.
Upon selecting the Query operation, ArcGIS Server will allow you to dynamically create/test and run custom queries, and provide you with variable output formats (HTML, json, KMZ, GeoJSON, and PBF). Due to GeoJSON’s flexibility, and ease of use in R, we will focus on GeoJSON query output formats throughout this exercise.
How do you know what fields/criteria to query on?
There are a couple of easy ways to determine criteria upon which to build your queries around. The first, and most immediate way is to simply examine the layer information at the REST endpoint which will show all field names and aliases. Below is an example of the fields present in the TMDL Watersheds layer:
The second way to really view the underlying data, is to leverage our web GIS apps to poke into the attribute and spatial data itself. Below are some examples of using the GIS Staff App to view and filter on specific attributes.
We will start by opening the attribute table for the TMDL Watersheds layer in the DEQ Data Layers list:
With the Attribute Table widget open, we can now build and apply validation-based filters to view the actual values in any field:
To view the unique values found in any field, simply select “Unique” in the “set input type” drop down. After a second, you will be able to view all of the valid values in whichever field has been specified. Here, we are looking at the Pollutant Name (POL_NAME) field, and viewing valid categories for that field.
Okay…let’s build a simple query and see it in R…
For this first example, we will query the DEQ REST service underlying the GIS Staff App. The REST url for this service is: https://gis.deq.virginia.gov/arcgis/rest/services/staff/DEQInternalDataViewer/MapServer
We’ll start by querying the WQM Stations layer (layer/ID = 104). We can see there is a wide range of fields that we can query this data layer on. For this first example, we will simply query for a single Station, with a Station ID (or WQM_STA_ID) of 2-FLA028.98.
So what’s going on below?…
We parse the URL to create a list where we can add all the parameters of our query - where, outFields, returnGeometry and f - with their appropriate values. As soon as the list is populated we use the function build_url() to create a properly encoded request. This request is the passed to the function st_read() to populate Stream_Stations. This Stream_Stations becomes both a simple features object and a data frame, i.e. it does contain both the geometry and the attribute values of the Station (in this case, a singular station). The interactive inlmisc::CreateWebMap() has various basemaps, a ID/Popup (on hover) based on the Station ID (can be changed to other fields), as well as the ability to zoom in/out:
#Read in GIS REST Data in geojson format
library(tidyverse)
library(httr)
library(sf)
library(leaflet)
library(inlmisc)
url <- parse_url("https://gis.deq.virginia.gov/arcgis/rest/services")
url$path <- paste(url$path, "staff/DEQInternalDataViewer/MapServer/104/query", sep = "/")
url$query <- list(where = "WQM_STA_ID = '2-FLA028.98'",
outFields = "*",
returnGeometry = "true",
f = "geojson")
request <- build_url(url)
Stream_Station <- st_read(request)
CreateWebMap(maps = c("Topo","Imagery","Hydrography"), collapsed = TRUE) %>%
addCircleMarkers(data = Stream_Stations,
color = 'red',
label = ~WQM_STA_ID )
As this data is now available for viewing, we can easily view/summarize/manipulate it further via summary, etc.
Querying multiple records:
Let’s look at the capability and performance of querying a series of records, based off of a more complex where statement. In this example, we are searching for all stations records attributed to the Appomattox River":
#Read in GIS REST Data in geojson format
url <- parse_url("https://gis.deq.virginia.gov/arcgis/rest/services")
url$path <- paste(url$path, "staff/DEQInternalDataViewer/MapServer/104/query", sep = "/")
url$query <- list(where = "WQM_STA_STREAM_NAME = 'Appomattox River'",
outFields = "*",
returnGeometry = "true",
f = "geojson")
request <- build_url(url)
Stream_Stations <- st_read(request)
CreateWebMap(maps = c("Topo","Imagery","Hydrography"), collapsed = TRUE) %>%
addCircleMarkers(data = Stream_Stations,
color = 'red',
label = ~WQM_STA_ID )
This works on lines and polygons as well!:
We can apply this same logic and format to query and spatially display polygon based layer. For example, we will query the TMDL Watersheds layer (Layer ID 72), searching for all Benthic Impairment TMDLs. We can also convert this data and read it into a data frame for further interpretation:
url <- parse_url("https://gis.deq.virginia.gov/arcgis/rest/services")
url$path <- paste(url$path, "staff/DEQInternalDataViewer/MapServer/72/query", sep = "/")
url$query <- list(where = "IMP_NAME = 'Benthic-Macroinvertebrate Bioassessments (Streams)'",
outFields = "*",
returnGeometry = "true",
f = "geojson")
request <- build_url(url)
Benthic_TMDLs <- st_read(request)
pal <- colorFactor(
palette = "Set2",
domain = unique(Benthic_TMDLs$WSHD_ID))
CreateWebMap(maps = c("Topo","Imagery","Hydrography"), collapsed = TRUE) %>%
addPolygons(data = Benthic_TMDLs,
color = ~pal(WSHD_ID),
fillColor = ~pal(WSHD_ID),
fillOpacity = 0.6,
stroke=5,
label = ~WSHD_ID )
Building more advanced queries
You can also add as many criteria as the input data allows. In this example, we are querying TMDL records with Sediment pollutants in the VRO region. When classified by status, we can quickly see that all qualifying TMDL records are approved.
url <- parse_url("https://gis.deq.virginia.gov/arcgis/rest/services")
url$path <- paste(url$path, "staff/DEQInternalDataViewer/MapServer/72/query", sep = "/")
url$query <- list(where = "POL_NAME = 'Sediment' AND REGION = 'VRO'",
outFields = "*",
returnGeometry = "true",
f = "geojson")
request <- build_url(url)
VRO_Sediment_TMDLs <- st_read(request)
pal <- colorFactor(
palette = "Set2",
domain = unique(VRO_Sediment_TMDLs$PROJECT_STATUS))
CreateWebMap(maps = c("Topo","Imagery","Hydrography"), collapsed = TRUE) %>%
addPolygons(data = VRO_Sediment_TMDLs,
color = ~pal(WSHD_ID),
fillColor = ~pal(PROJECT_STATUS),
fillOpacity = 0.6,
stroke=5,
label = ~WSHD_ID ) %>%
addLegend(data = VRO_Sediment_TMDLs,
pal = pal,
values = ~PROJECT_STATUS,
title = "Legend",
position = 'topright')
Armed with your new confidence and familiarity with querying and interacting with ArcGIS Server web services… go out into the world and see what kinds of data you can query, summarize and analyze!
